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ABSTRACT 

The kinetic SZ effect, which is the dominant CMB source at arc-minute scales and 
v ~ 217 Ghz, probes the ionized gas peculiar momentum up to the epoch of reion- 
ization and is a sensitive measure of the reionization history. We ran high resolution 
self-similar and ACDM hydro simulations and built an analytical model to study this 
effect. Our model reproduces the ACDM simulation results to several percent accu- 
racy, passes various tests against self-similar simulations, and shows a wider range of 
applicability than previous analytical models. Our model in its continuous version is 
free of simulation limitations such as finite simulation box and finite resolution and 
allows an accurate prediction of the kinetic SZ power spectrum C/ . For the WMAP cos- 
mology, we find Z 2 C//(2tt) ~ 0.91 x 1(T 12 [(1 + z roion )/10] a34 (V5000)°- 23 -°- 015 ( z — " 9 > 
for the reionization redshift 6 < z rcion < 20 and 3000 < I < 9000. The corresponding- 
temperature fluctuation is several /iK at these ranges. The dependence of C; on the 
reionization history allows an accurate measurement of the reionization epoch. For the 
Atacama cosmology telescope experiment, C; can be measured with ~ 1% accuracy. C/ 
scales as (f2(,ft,) 2 cr|~ 6 . Given cosmological parameters, ACT would be able to constrain 
Zreion with several percent accuracy. Some multi-reionization scenarios degenerate in 
the primary CMB temperature and TE measurement can be distinguished with ~ 10<r 
confidence. 

Key words: Cosmic microwave background-theory-simulation: large scale structure, 
intergalactic medium, reionization, cosmology 



1 INTRODUCTION 

Reionization is a critical phase in the evolution of the intergalactic medium (IGM). Recent constraints on the reionization 
epoch based on Lyman-a absorption and CMB appear to be in conflict. The WMAP measurement of the CMB temperature 
and polarization fluctuation constrains the Thomson optical depth to the last scattering surface r ~ 0.17 ± 0.04. This large 
optical depth sugg ests an early reionization at z ~ 17 if the universe is re-ionized suddenly and remains completely ionized 
after that l)Kogut et alib oOS: ISoergel et all2o'o^) . This inference of the reionization redshift and optical depth are very robust: 
Thomson scattering is polarization dependent, and converts a temperature quadrupole into a linear polarization. The observed 
~ 2\xK temperature polarization cross correlation are just the product of the optical depths and the local CMB quadrupole, 
Qrms.ps = 20/iK, which is known to high precision. Given the optical depth, and the known baryon content of the universe, 
one obtains a minimum redshift to accumulate the sufficient total column density of electrons. 

But, detection of a Gunn-Peterson absorption trough in a SDSS quasar at redshift z = 6.28 sets a constraint on the 
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rcionization redshift at z ~ 6 teecker et alJbOQlh . The quasar radiation is seen to be fully absorbed at a Lyman-a redshift 
larger than 6.1. This corresponds to a lower bound of neutral hydrogen column density, corresponding to a neutral fraction 
greater than a percent. At face value, it would appear to be possible to have a small neutral fraction while the majority of 
the electrons are ionized. But if ionization occurs through the UV radiation from discrete sources, this radiation must be able 
to penetrate the IGM to ionize it. It is generally not possible for the radiation to pass through an opaque medium to ionize 
material behind it. One would generically expect the universe to be either mostly neutral, or fully transparent with neutral 
fractions less than 10 -3 . Furthermore, in order to keep the IGM temperatu re T ~ 2 x 10 4 K at z ~ 2-4 required by the 
Lyman-a forest observation, a late time reionization must occur after z = 10 faui fc Haimanll2003l) . 

Exotic possibilities exist to reconcile a sudden reionization with both WMAP polarization and SDSS quasar data. The 
ionization could arise from decaying dark matter, which would be uniformly distributed. Or it could have ionized multiple 
times, and recombined in-between. All these arguments lead to more complicated reionization scenarios. 

To make observational progress on this conundrum, more information about the fraction and dynamics of free electrons at 
z ~ 10 is needed. The primary CMB temperature measurement depends on the reionization history only through r and thus 
has an intrinsic deg eneracy. A precision CMB E-mode polarization measurement is in principle able to break this degeneracy 
JHu fc Holderll2 003V But the weak polarization signal makes such measurement unfeasible in the near future. Due to the huge 
Lyman-a optical depth which must be present in any multi-reionization scenarios, the Lyman-a forest can only detect the last 
reionization a nd makes it an indirect probe of the reionization history. The 21 cm e mission and absorption backgrounds in 
the dark age l|Tozzi et alJ boOol llliev et a ll bool Eool Ichen fc Miralda-EscudeH 200 3) would constrain the upper limit of the 
first reionization redshift. The contamination by free-free emission of ionizing halos joh fc MacklbOQ.'t ) may be distinguished 
by their spectral features. Such a measurement may be possible with a square kilometer array, but its construction is still far 
in the future. 

Free electrons resulting from reionization also scatter off CMB photons through Compton scattering and result in the 
kinetic Sunyaev Zel'dovich (SZ) effect. Compton scattering keeps information from high redshifts because it does not depend 
on redshift and is not affected by distance or the expansion of the universe. Since free electrons at high redshifts have higher 
number density, their contribution to the kinetic SZ effect is larger. So the kine tic SZ effect is a sensitive measure of the 
reionization process at high redshifts. Next generation CMB experiments such as lthe Atacama Cosmology Telescope! l|2003^ 
are making its precision measurement feasible in the foreseeable future. Once we have an accurate model of this effect, the 
reionization history will be precisely constrained. It is the goal of this paper to build such a model and to quantify the power 
of the kinetic SZ effect to constrain the reionization history. 

Unlike most radiative processes, the physics of the kinetic SZ effect is very simple. In contrast to models of the Lyman-a 
forest, the ionized gas peculiar momentum does not require an accurate understanding of the gas state such as metalicity, 
temp erature and ionization equilibrium. In the linear regime, the exact model is known as the Vishniac effect ijVishniad 
1987). Th e non-linear regime is intrinsica lly more difficult and theoretical works in th e literature have yet to converge quan- 
titatively faJbood Ida Silva et aJboOlallbt iGnedin fcJaffelboOll : ISpringel et allboOll : IValageas et"alll200lt iMa fc FrvlEool 
IZhang. Pen fc Wanell2002l) . 

Hydrodynamical simulations are best equipped to capture the nonlinear IGM physics. Large box sizes are required to 
probe large-scale power while high resolution is needed to resolve small-scale, nonlinear structures. Most simulations to date 
have sacrificed one for the other. But the kinetic SZ effect requires the capturing of both large and small scale structures 
and puts a strong requirement on the simulation power. Limited computational resources introduce numerical artifacts which 
must be quantified. In contrast, analytical models can be constructed to span all relevant scales continuously, but they are 
often ad hoc procedures and must be tested and calibrated against simulations. In this paper, we run Eulerian hydrodynamical 
simulations with box sizes of 50 and 100 Mpc//i to model the evolution of the IGM and measure the gas momentum for various 
redshifts. Numerical limitations such as finite simulation box size, finite resolution, and non-linearity are studied by running 
different resolution simulations, ranging from 256 3 up to 1024 3 fluid elements. Self-similar simulations are also used to quantify 
the numerical limitations. We then build an analytical model covering both linear and nonlinear regimes, motivated by the 
Vishniac effect and simulations. Our model is free of numerical artifacts and captures the main effect of the non-linearity at 
the same time. 

In J5|we review the formalism for the KSZ effect. In S^we describe the numerical simulations, including their results and 
their limitations. Motivated by the gas hydrodynamics and simulations, we build an analytical model in This model is 
tested against simulations in In f||3we calculate the KSZ power spectrum using our analytical model and in 33 we discuss 
the potential of the KSZ effect to constrain the reionization history. 



2 THE KINETIC SZ EFFECT: FORMALISM 



The kinetic SZ effect directly probes the ionized ele ctron peculiar momentum p = (l + <5 e )v through Compton scattering. The 
resulting CMB temperature fluctuation is given bv lSunvaev fc Zel'dovichl Jl98d) : 
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Q(") = = Xen e a T ^—^exp[-T(z)]adx= ^—^exp[-r(z)]dT. (I) 

-1CMB J C J C 

Here, n e is the mean electron number density and \ e is the ionization fraction, or is the Thomson cross section, h is the 
direction on the sky and x is the comoving distance. r(z) — f^ x £ n £ aTaAx is the Thompson optical depth, p can be decomposed 
into a curl- free or gradient part p_e satisfying V x pg = and a divergence-free or curl part p_e satisfying V • ps =0. The 
gradient term cancels out when integrating along the line of sight and has no contribution to the kinetic SZ effect. Thus the 
only part which contributes to the integral Q is ps. We define the correlation functions £b,e(t) = (ps,s(x) ■ pb,e(x + r)) 
and the corresponding power spectra P| = {|p_B.E(k)| 2 }. Throughout this paper, we use A%{k,z) = fc 3 Pl(fc)/(27r 2 ) as 
the gas momentum curl part power spectrum instead of p%(k). A 2 B {k, z) is a detailed description of the IGM state, whose 2D 
projection is the kinetic SZ power spectrum Ci (Eq.[3Jl. 

Adopting Limber's approximation and the fact that (|p s (k) -n\ 2 ) — ^Pg(k), we obtain the kinetic SZ angular correlation 
function: 



w{6) ~ (^J cose J dx[an eX e] 2 cxp[-2r(z)] J ^(y/x^+y^dy (2) 

f ^ T ' 



y — j / dx[an e x e ] 2 exp[-2r( Z )] / ^/xW + y*)dy. 



Here, :r ro is the comoving distance to the reionizati on epoch. We have s et \e = Xe, n amely, omitted the p atchy reionization 
effect 1 following the suggestion of both simulations jGnedin &Jaffell200lt) and theory JValaeeas et aljfeooil) . T he last approx - 
imation introduces an error less than 0.2% for 6 < 1°. Thus, the usual expression of the Limber's equation iPeacockHl999l) 
still holds: 

Ci = (2^)3 [^7^] 2 >< j ( ZrC (l + 2 ) 4 X^A|(fc,,)U. i/ ,exp[-2r W ] a;W ^^. (3) 

ps is zero in first order linear perturbation theory. Thus, the kinetic SZ effect is an intrinsically non-linear process which 
must be investigated by either hydro si mulations or n on-linear theories such as high order perturbation theory. After the 
pioneering work of the Vishniac ef fect lYi^hniacJIl^sj), w hich is a robust pred i ction of the kinetic SZ effect in the linear 
regi me, various authors Eul 120001: Ida Silva et, alJ foffllallbl lOnedin fc.Taffel l200lt ISoringel et"alll200lt IValageas etafl l200ll : 
iMa fc Fnl200i Izhang. Pen fc Wanj |2002l have studied the kinetic SZ effect in the non-linear regime with analytical models 
and simulations, but there are still some unclear key issues. Three key issues we will address in this paper are (1) How well 
do simulations capture A% and thus the large scale power spectrum of the kinetic SZ effect at large and small scales? 
has contributions from both large and small scales. Simulations have cutoffs at both scales due to the finite simulation box 
size and resolution, respectively. These cutoffs affect simulation results and need to be quantified. Then, how to quantify 
and overcome these effects to produce reliable results? (2) What is the effect of the non-linearity to A%? How to model this 
effect? (3) How observable is the effect, and to what accuracy can it be used to model measured reionization by upcoming 
experiments? We will address these issues using simulations and our analytical model f tHp . 



3 HYDRODYNAMIC SIMULATIONS 

We ran cosmological hydrodynamical simulations using a new Eulerian cosmological hydro code <Trac fc Penll2003allrl . This 
Eulerian code (hereafter TP) is based on the finite- volume, flux-conservative total variation diminishing (TVD) scheme that 
provides high-order accuracy and high-resolution capturing of shocks. The hydrodynamics of the gas is simulated by solving 
the Euler system of conservation equations for mass, momentum, and energy on a fixed Cartesian grid. The gravita tional 
evolution of the dark matter is simulated using a cloud- in-cell particle- mesh (PM) scheme llHocknev fc EastwoodllT988|) . 

The robustness of the TP code has been tested b y comparing the ev olution of the dark matter and gas density power 
spectra from the simulations with the fitting for mula oflsmith et al. I J2OO2I) . We also performed a code comparison by running 
the same initial conditions using the MMH code ijPenllf 99SF) . which combines the shock capturing abilities of Eulerian schemes 
with the high dynamic range in density achieved by Lagrangian schemes. Power spectra are computed using FFTs. We find 
good agreement at all relevant scales and redshifts for both comparisons. 

The KSZ effect is a non-local problem and has contributions from both large-scale and small-scale power. This makes it a 
challenging problem to model numerically because high resolution is needed at all scales. Eulerian schemes provide high mass 
resolution at all scales, unlike Lagrangian schemes based on SPH which achieve high dynamic range in high-density regions 
but do poorly in low-density regions. Eulerian schemes are ideal for simulating the evolution of the IGM to model the thermal 
and kinetic SZ effects and the Lyman alpha forest, because of their high speed, superior mass resolution, shock-capturing 
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abilities. Furthermore, Eulerian algorithms are computationally very fast and memory friendly, allowing one to optimally use 
available computational resources. 

We ran a t otal of six simulation s with the best fit WMAP-alone cosmology fl m = 0.268, S1a = 0.752, Qt = 0.044, h = 0.71, 
and as = 0.84 JSpergel et al.ll2003h . Two box sizes of 50 and 100 Mpc/fa are used and for each we perform the simulations 
on fixed grids with 256 3 , 512 3 , and 1024 3 cells. The ratio of dark matter particles to fluid elements is 1:8 in all six ACDM 
simulations. In the highest resolution simulation, we achieve a spatial resolution of Ax ~ 50 kpc//i and a dark matter particle 
mass resolution of Am ~ 7 x 10 7 Afp. The initial con ditions are generated by sampling from an initial power spectrum 
computed using CMBFAST dSeliak fc Zaldarriagairi99d) . The simulations are started at a redshift of z = 100 and evolved 
down to z = 1, with data outputs at z = 20, 10, 6, 3, and 1. The KSZ effect, for angular scales I > 1000, has contributions from 
electron scattering primarily from the z > 1 IGM and therefore, we do not push the simulations below z — 1 for the purposes 
of this paper. Furthermore, the final prediction of the kinetic SZ effect is calculated using our analytical model, instead of 
using these simulations directly, where the contribution from z < 1 IGM is included. The highest resolution simulation takes 
approximately 700 time-steps to evolve from z — 100 down to z = 1. On a GS320 Compaq Alpha server with 32 cpus and 
total theoretical peak speed of 32 Gflops, the run takes approximately 2 days. 

The peculiar gas momentum is decomposed into the rotational (B) and irrotational terms (E) in Fourier space. In Fig.0 
we plot the simulated dimensionless As(fc, z)k/H(z) at z = 1 and 10 and compare them to the linear Vishniac effect, which 
is a robust prediction in the linear regime (J1J. H(z) is the Hubble parameter. As(fc, z)k/H(z) is the velocity of a wave 
expressed in units of wave lengths per expansion time. It reaches ~ 1 when the overdensity reaches the mildly nonlinear value 
5 ~ 1. We note that our convention differs from the literature that we plot the root-mean-square standard deviation instead 
of a variance. 

Three key conclusions can be drawn from our simulations. (1) Finite simulation box size causes A%(k, z) to lose power at 
large scales. At redshift z — 10, where linear theory is a good approximation, we would expect the simulation results to agree 
well with the Vishniac prediction. But it turns out that they agree only at a very limited k range. More counter-intuitively, 
the discrepancy diverges toward large scales. The simulated A%(k,z) loses power toward large scales. 50 Mpc//i simulations 
lose more power at large scales than 100 Mpc/ft simulations at the same scales. This suggests the effect of the finite simulation 
box size. Simulations cut off density fluctuation modes larger than half simulation box size and thus affect A% through its 
non-local dependence on the density and velocity fields. (2) Finite simulation resolution causes A%(k, z) to lose power at small 
scales (less than ~ 10-20 grid separation). A%(k,z) truncates at small scales and the truncating scale fct run c keeps increasing 
with respect to the simulation resolution. Gas pressure is expected to decrease power at very small scales. If this truncation is 
caused by the gas pressure, we would expect fct run c to decrease with the simulation resolution since higher resolution results in 
higher gas pressure. Since we observe opposite behavior, we believe that this truncation is caused by the simulation resolution, 
which predicts that fc trU nc to be roughly proportional to the resolution, as apparent in figure. (3) The non-linearity increases 
A%(k, z) with respect to the linear prediction. This is visible in the top panel of figure. Q 

We further investigate these issues by self similar simulations. We ran two self similar 512 3 simulations with n = — 1 and 
n — —2, respectively. Q m = 1.0. fif, = 0.164 is adopted such that fij,/fi m is the same as the WMAP cosmology . The effects of 
finite simulation box size and finite simulation resolution scale as expected. The non-linear effect turns out to be more subtle. 
For n = — 1, simulated A% is smaller than the Vishniac prediction in the non-linear regime. The non-linearity decreases A|. 
For n = —2, simulated A% is larger than the Vishniac prediction. Thus the non-linearity increases A%. To quantify simulation 
artifacts and the non-linearity effect, we build an analytical model in Therefore we postpone the discussion and explanation 
of the self similar simulations until in f0 



4 THE ANALYTICAL MODEL 

In this section, we will extend the Vishniac result to the grid space to model the simulation effect ( jHfl and to the non-linear 
regime to model the non-linearity effect fi |4,2t . 

4.1 Linear regime 

In the linear regime, the motion of IGM is determined by the gravitational potential <j>. Thus v o c V4> and i s curl -free. Since 
V 2 <f> oc S, v(k) oc 5 k k/k 2 . From this relation, A^ in the linear regime was originally calculated bv lVishniad <ll987l) . With the 
correction of the gas window function W g , we obtain 

Ai(fc) = ^\ {<1 ^ 2 J P L (k 2 ,z)P L (\k - k 2 \,z) [W g (\k - k 2 |)/3(k,k 2 ) + W g (k 2 )f3(k,k- k 2 )] 2 d 3 fc 2 . (4) 

D is the linear growth factor of the dark matter over-density, D is its time derivative and a is the scale factor of the 
universe. Ph{k, z) is the linear dark matter power spectrum. Gas pressure originated from various non-gravitational heating 
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Figure 1. The simulated Aj(i;)fc/ H(z) at Z = 1 and 2 = 10 for four of our simulations. We are plotting the standard deviation As and 
not the variance A B which is usually plotted. Comparing to the Vishniac prediction, which is accurate at linear scales, simulations lose 
power at these scales. The 50 Mpc/h simulation loses more power than the 100 Mpc/h simulations and suggests the effect of the finite 
simulation box size. At small scales, the simulation resolution effect is shown in the large k behavior of the three 100 Mpc/h simulations. 
Simulated As(fc) begins to lose power at less than 10-20 grid separation. Ag(fc)fc/ H{z) ~ 1 signals the linear-nonlinear transition region. 
Non-linear effects increase A B (k) relative to Vishniac theory for the WMAP cosmology we adopted. We also notice that the sign of the 
correction depends on the slope of the power spectrum. 



processes smooths out the gas density at small scales. This effect can be modeled b y a window fun ction i y Q (fc) su ch that 
PJk) = Wh k )PpM (k), which can be choose as a Gaussian window function (see e.g. Irlui fc Gnedlnl Jl997t ) : iGnedinl J2OO0I) : 
iGnedin et alJ teoOotO . In an adiabatically evolving universe, gas traces dark matter down to very small scales and W g would 
be effectively unity. One can model the simulation resolution by introducing an effective grid smoothing window W g . This 
procedure works well but the modeling of W g is arbitrary. For simplicity, we omit this gas window effect. Since the gas 
peculiar velocity v g is mainly determined by the dark matter gravitational potential field, we have assumed that v g traces the 
dark matter peculiar velocity at all scales. The kernel /3(k, k 2 ) = (k 2 — j^! )/k|. This form originates from the gradient 
operators in v a V<^> and S7 2 cj> oc S. Eq. @ explicitly preserves the symmetry between ki = k — k 2 and k2. 

Eq. (J1J is exact in the linear regime and an ideal simulation should be able to reproduce these results. But we notice from 
this equation that A% has contributions from both large and small scales. For a real simulation, the simulation box size and 
resolution cut off contributions from these scales and the simulation result may deviate from Eq. IgJ. To address this issue 
and to estimate the limitation of simulations, we calculate A% on discrete grids. For a simulation with box size L, resolution 
N and the periodic boundary condition, the Fourier component of the gradient operator V is V = — i sin(kAL)/AL = 
— i {sm(k x AL)/AL, sin(k y AL)/AL, sin(fc z AL)/AL}. Here, AL = L/N. Comparing to V = — ik in the continuous case, we 
obtain the grid version of A 2 B : 

A|(fc) = ^\(ab/DfL z J2 «)A(|(k - k 2 ) P |, z) [W„(|(k - k 2 ) P |)/fc(k,k 2 ) + W g (k 2 )f3(k, (k - k 2 ) P )] 2 . (5) 

k 2 

Here, /3 G (k,k 2 ) = /3(k -> sin(kAL)/AL, k 2 -> sin(k 2 AL)/AL), k x , y , z ,k 2x ,y : z = 2n/L x {-N/2, -N/2 - l,...,N/2 - 1). 
(k — k 2 ) P means k — k 2 under the periodic boundary condition. Namely, k — k 2 needs to be converted to the same range 
as k 2 by the periodic boundary condition. This grid version represents simulation discretization effect at the first order 
approximation. 
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4.2 Non-linear regime 

In the non- linear regime, the perturbation theory does not hold and Eq. @J & © break down. Irlul J2OO0I) : iMa fc Frvl (120021) 
argue that the velocity field is less non-linear than the density field and one may substitute Pl(/s, z) in Eq. Q introduced by 
the density field with the corresponding non-linear density power spectrum PMh{k,z) and keep the other Pl introduced by 
the velocity field unchanged. But the generation of curl in the velocity field by shell crossing in the non-linear regime weakens 
this argument. In the presence of curl, V x ps = V<5 x v_b + V<5 xvs + (l + 5)Vxvs. The first term i n the right side of this 
equat ion is well described by the linear Vishniac effect and may be well described by the method of Irlul J200(J) : livia fc Frvl 
( 2002) in the non-linear regime. But the remaining terms bring non-negligible contribution in the non-linear regimes. In 
our simulations, we find that in the highly non-linear regime, velocity field reaches equi-partition an d A ? , = 2 A 2 rrj . This 
sugges ts that the contribution from vs may be dominant in the non-linear regime. Thus, the method of lHul l)2000h : lMa fc Frvl 
( 2002) underestimates A| . The generation of v_b through shell crossing involves complicated gas hydrodynamics and makes 
the analytical understanding of v_b intractable. Its statistics in principle can only be probed by simulations. But since the 
generation of vs is caused by the non-linearity, to model its effect, one may phenomenologically substitute the linear density 
power spectrum in Eq. Q introduced by the velocity field with the corresponding non-linear density power spectrum. This 
procedure at least qualitatively captures the effect of vb generation. Thus we obtain the final continuous form of A^(fc, z): 

A|(fc, z) = J^l {a ^/ J P NL (k 2 ,z)P NL (\k - ka|,2) [W g (\k - k 2 |)/3(k,k 2 ) + W 9 (fe,)/9(k,k- k 2 )] 2 d 3 fc 2 (6) 
and the corresponding grid form 

A|(fc) = ^±(ai)/D) 2 L 3 PN L (k 2 , z)P NL (\(k - kaU*) [W g (\(k - k 2 ) p |)/3 G (k,k 2 ) + +W fl (fa)/J(k, (k - k 2 ) p )] 2 . (7) 

k 2 

In our mod el, the nonlinearity e ffect is explicitly quadratic and larger than that of lHul (200n^: lMa fc Frvl |2002). We will adopt 
the code of ISmith et al. I ^2002) to calculate PMh{k,z). 

In 50 we will test our model against simulations. Our logic is, if the grid version of our model, which captures key 
simulation artifacts, gives a good description of the simulation, then, its asymptotic case, namely the continuous version of 
our model would represent an ideal simulation, be free of these simulation artifacts and describe the real A%. 
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Figure 3. The comparison between our model and the simulations with WMAP cosmology. The vertical axis is the dimensionlcss 
A B k/ H(z), instead of the usually plotted. Two sets of data points are our 1024 3 , 100 Mpc/h (open square) and 1024 3 , 50 Mpc/h 
(filled hexagon) results. For clarity, we only show a few data points. As comparison, we show grid version predictions (dash lines), 
continuous version predictions (solid lines) of our model and the the Vishniac predictions (dot lines). Our grid predictions agree with 
simulations at all redshifts and all scales that simulations can reliably probe. 



5 COMPARISON BETWEEN ANALYTICAL MODEL AND SIMULATIONS 

We first compare the grid version prediction of our model Eq. J7J with simulations. Self similar simulations are idea l to test 
the ef fects of box size, resolution and non-linearity due to the scaling relation between redshifts and scales (See e.g. IPeeblesI 
lllflSftl V The grid prediction of our model fits simulations very well at all reliable scales. The cut off due to the finite simulation 
box size can be seen from the large scale behavior of A% and the resolution effect can be seen from the small scale behavior 
of A| (Fig. HI. The depend ence of the non-linearity effect on the initial density power index n is naturally explained by our 
model. The effect of the non-linearity to the density power spectrum depends on n. For n = — 1, the non-linearity decreases the 
non-linear density power spectrum with respect to the linear density power spectrum. Our model then predicts a decreasing 
in A| with respect to the Vishniac prediction by the non-linearity. For n — —2, it is the opposite case. For the WMAP 
cosmology, the effective power index in simulation scales is around —2 and an amplification of A B by the non-linearity is 
expected, as shown in Fig. 

Since our ultimate goal is to provide an accurate description of A B for our real universe, we test our model in detail 
against the 1024 3 simulations with WMAP cosmology. In order to do that, we need to estimate the variance of A B (k). In the 
linear regime S is Gaussian distributed and ps(k) cx S5^, we obtain {p%(k)} — ^(ps(k)) 2 and estimate the uncertainty of 
the simulated A 2 B by this relation. Our result (Fig. shows that the grid calculation agrees with simulated A 2 B at all output 
redshifts and the full k range that simulations can reliably probe. Our model passed the test of all self similar and WMAP 
simulations we have run and shows a wider range of applicability than all previous analytical models. Thus we believe that 
the grid version of our model is a suitable tool to describe simulated A 2 B . Its continuous version (Eq.|SJ then should be able 
to describe an ideal simulation with infinite box size and resolution and produce the real A%. 

Comparing our analytical prediction and simulations, we concluded that the simulated A| loses power at large scales, 
gains power at small scales due to non-linearity until approaching the resolution limit. Thus it has a steeper slope. We address 
several simulation limitations here. (1) The finite simulation box size causes the simulated A 2 B to lose power at large scales. 
From Eq. 0, A| has contributions from all scales. The finite box size cuts off the density fluctuation modes with size larger 
than half box size. It is these missing modes causing the simulated A 2 B to lose power at large scales. The smaller the box size, 
the more A 2 B power is lost, as can be seen from Fig. [3] (2) The finite resolution causes A 2 B to lose power at less than about 
10-20 grid separation. (3) The periodic boundary condition causes the simulated A B (k) to gain power. This can be seen by 
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0.1 1 10 100 

A 2 (k) 

Figure 4. The correlation coefficient r between <5 and 5. r = A 2 ^/ ^A|A|. In the highly non-linear regime, 6 and S are barely correlated 
and r 2 roughly scales as 1/A|. The result shown is for our WMAP simulation with 50 Mpc/h box size and 1024 3 grid elements. 



the following argument. For the cosmology we adopt, P(k) is a decreasing function for k > 0.02 Mpc/h. Comparing Eq. £j) 
and Eq. Q, the periodic boundary condition produces a smaller k — k2 and thus a larger P(|k — k.2 1) for those k — k2 which 
are beyond the simulation k range. The case is the same for the gas window function and is similar for the kernel function 
/3(k,k 2 ). But this effect is negligible comparing to the first two. 

Our modeling of the non-linearity effect works surprisingly well. For WMAP co smo l ogy, t h e non-linearity do es increase 
A 2 B significantly (fig. |3J. As expected, the non-linearity effect is larger than what iHul 1120001 ) : iMa fc Frvl 1120021) predicted 
since their models omit the contribution of v_g. Our model passes the test of all of our simulations and suggests that it not 
only captures the contribution of v_b qualitatively but also quantitatively. Though it is hard to explain why it works well 
quantitatively, we believe that it gives a reasonable description of A 2 B in the non-linear regime. 

Our model could be understood from the equi-partition of p in the non-linear regime, which states A 2 B = 2A%, as found 
in our simulations, and the mass conservation equation, which states S + V-p^/a = 0. Defining the cross correlation coefficient 
r = A^/y^A 2 A? between 5 and 5, one then has 

A 2 la 2 A 2 \k) 2 

(8) 



2 k 2 r 2 A 2 (k) 



r in the linear theory and the stable clustering regime can be calculated analyticall y. But since possibl y due to the ongoing 
merger process, the stable clustering may not hold even in highly non-linear regime dSmith et al. 1 ^002 V for most non-linear 
regimes, r is hard to estimate, even in an interpolation way. In our simulations, we find that in highly non-linear regime, 5 
and 5 becomes weakly correlated . r 2 roughly scales as 1/A 2 (Fig. QJ. This scaling gives A% oc A 4 , which is consistent with 
our previous model (Eq.|Hland|7J. 

Simulation limitations put strong requirements on the simulation power. For example, from Eq. 10, we estimate that, 
in order to produce a reliable A 2 B at 0.1 h/Mpc < k < 10 h/Mpc, a hydro simulation with a box size L ^ 400 Mpc/h and 
resolution > 4096 is required, which is not realistic at present. But the continuous version of our analytical model (Eq. 

is free of these simulation limitations and is able to produce a reliable A% at all k and z range relevant to the kinetic SZ 
effect. Hereafter, we will adopt Eq. © to calculate the kinetic SZ power spectrum. 
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Figure 5. The kine tic SZ Ci dependence on scales and reionization redshift. The primary CMB is the lensed CMB calcu- 
lated by CMBFAST fceliak fc Zaldarriagal [l996l) . For 6 < z reion < 20 and 3000 < I < 9000, l 2 Ci/(2n) ~ 0.91 X 10" 12 [(1 + 
^cion)/10] - 34 a/5000) - 23 -°- 015 ( 2 -o„-9). 



6 THE KINETIC SUNYAEV ZEL'DOVICH EFFECT 

The kinetic SZ power spectrum Ci is very sensitive to the reionization history. To see this, we estimate its dependence on the 
reionization redshift. For simplicity, we assume that the universe is reionized at once at z rc ion and remains completely ionized 
after that. Since r C 1, we omit its effect temporarily for the estimation. First, we consider the linear case. Combining Eq. 
@ & @, we obtain 

d(z Icion ) oc I rC '°° (1 + z) 4 (aDD) 2 ( X (z)) 1 - nB ^dz oc J IO1 ° a (^-D) 2 ^/(n (l + z)3+n v )( X (z)) 1 - nB dz (9) 

"(1 + z) (1 -" B)a -°' 5 d 2 oc (1 + Zrcio n) (1 - nB)Q+0 - 5 for z > 1 

Here, ns is the power index of A%(k, z) at k = l/\(z). For / ^ 10 3 and z > 1, ns < 1. We approximate xi z ) oc (1 + z) a . For 
z ^ 2, < a < 0.5. For z ^ 10, < a < 0.1. When z > 2, D oc a and ^ oc 1. Then, for z > 2, Ci diverges toward high z rc ion 
with a scaling relation Ci oc (1 + z rc i on )~ ' 5 . The non-linearity increases the contribution from low redshifts and makes the Ci 
redshift dependence weaker. The damping caused by the optical depth r causes further suppression. We show Ci calculated 
for various Zreion by Eq. © & © in Fig. Qj] The Ci behavior has three distinct regions. (1) I <C 1000. C; in this region is 
mainly contributed by low redshift linear regions. Its dependence on z rc i n is weak. (2) / 3> 10000. C; in this region is mainly 
contributed by low redshift highly non-linear regions. Its dependence on z rc ion is also weak. (3) 1000 < / < 10000. Ci in this 
region has large contributions from the high redshift universe and thus has the highest sensitivity to z rc io n . As we will discuss 
in J7| this region is also observationally accessible in the near future. For 6 < z rc i on < 20 and 3000 < I < 9000, our result can 
be fitted by 

l 2 Ci/ (2tv) ~ 0.91 x 10" 12 [(1 + 2 rcion )/10] ' 34 (Z/5000)°' 23 -°' 015(z — " 9) . (10) 

The error introduced by this fitting formula is less than 2%. The corresponding temperature fluctuation AT is then 

AT; ~ 2.60[(1 + z rcion )/10] - 17 (V5000)°' 115 ^ ' 0075(z ™ i °"- 9 VK. (11) 

For the WMAP favored z rc i on = 16.5, we show Ci predicted by our model and by the corresponding Vishniac effect in 
figED The non-linearity increases the power spectrum at small angular scales. At I ~ 4000, this amplification reaches 2. The 
non-linearity then produces a broad peak in the power spectrum, which extends from I ~ 3000 to I ~ 10 with an amplitude 




10 Pengjie Zhang, Ue-Li Pen & Hy Trac 



10 



-9 



10 



10 



^UO" 11 



10 



10 



10 



-12 



-13 



-14 



= III 

- 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 III- 

Primary CMB 








\ KSZ \ 


— 

// 
// 

= ~ / / 


\ VishrTi-aQ V 


— // 

- // 
// 

V 

1 III 


\lensed I 

Mill 1 



100 



1000 



10 4 



10 E 



Figure 6. The kinetic SZ effect power spectrum (solid line). The WMAP cosmology with tt m = 0.268, Ha = l-^o, &b = 0.04, h = 0.71, 
erg = 0.84 is adopted. We assume that the universe suddenly ionized at z = 16.5 and remain completely ionized after that. To show 
a compari son, we show the primary C MB power spectrum and the Thermal SZ power spectrum, which is adopted from a 512 3 MMH 
simulation Zha ng. Pen &; Wand 120021) and has been scaled to the WMAP cosmology by the scaling relation C; <x cr^ 7 . The dash line is 
the Vishniac prediction. The non-linearity increases the kinetic SZ power spectrum by a factor of 2 at I ~ 4000. 



AT ~ 2.7/iK. The large I behavior will be certainly affected by the gas window function W g (k). Its effect can be easily 
incorporated, as discussed in H4.ll Since it mainly affects the I > 10 4 region, where it is hard to observe, we will not discuss 
its effect here. 



7 THE POWER FOR THE KINETIC SZ EFFECT TO CONSTRAIN THE REIONIZATION HISTORY 

For Zrcion ~ 7, as favored by the Lyman-a absorption observation versus 2: rc i on ~ 17, as inferred from WMAP, Ci would differ 
by 45%. This sensitive dependence of Ci on the reionization history opens a window for future CMB observations to constrain 
the reionization history and break the degeneracy met in the primary CMB argument. The power of the kinetic SZ effect to 
constrain the reionization history depends on the accuracy of the kinetic SZ effect measurement. In this section, we take the 
Atacama cosmology telescope (ACT) as our target to address this issue. 

ACT will operate at three bands, 145 Ghz, 225 Ghz and 265 Ghz, measure over 100 square degree sky with noise error 
on — 2/iK per 1.7 x 1.7 pixel. At the 225 Ghz band, the thermal SZ signal is negligible and we take this band to estimate 
the sensitivity of ACT to measure the kinetic SZ effect. The noise power spectrum is given by Cn = Anf p i x a N . Here, / p i x is 
the fractional sky coverage of each pixel. The error in the kinetic SZ power spectrum measurement comes from three sources: 
instrumental noise, primary CMB and cosmic variance. If the primary CMB power spectrum and the noise spectrum are 
known and can be deducted, we have then 

Ad = I (y + 2) + 2(C N /C l /W N ) 2 + 2(C C mb/C 1 ) 2 

Ci V (2Z + 1)AZ/ Sky ■ 1 ' 

Wn(1) is the Fourier transform of the natural beam function. For the 225 Ghz band, the ACT resolution is 1.1 . We then 
approximate the beam function as a top-hat function with an effective radius y 1.1' x l.l'/w- We normalize it to have 
Wn(1 — ► 1) — > 1. /sky is the fractional sky coverage of the whole survey area, which we adopt as 100 square degree, r) is a 
measure of the non-Gaussianity, which will be taken as 105/9 — 3 = 26/3, as in the 3D case obtained from the second order 
perturbation theory ( H4.H . Due to the projection, the 2D kinetic SZ non-Gaussianity would be weaker and the error obtained 
above should serve as an upper limit. 
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Figure 7. The sensitivity of ACT to distinguish reionization scenarios. The three unlabeled lines are rcionization scenarios A, B and 
C, from bottom up. The error bars are the estimated ACT errors for the reionization scenario B. At I ~ 6000, the error is ~ 3%. B 
and C are degenerate in the primary CMB measurement but they can be distinguished by ACT with > 3<r confidence. One can bin the 
data points at 3000 < / < 8000 and do not lose useful information since the power spectrum is almost flat at this range, thus one can 
reach 1% accuracy in the C; measurement and B and C can be distinguished by lOcr confidence. For the single reionization scenario, the 
rcionization redshift can be constrained with better than 3% accuracy. 

We consider three reionization scenarios: (A) the universe reionized once at z rc ion = 7, (B) the universe reionized com- 
pletely once at z rc i on = 16.5, (C) the universe reionized completely at z = 21 and became neutral at z — 13 and reionized 
again at z = 7. The last two has roughly the same optical depth r = 0.174 to the last scattering surface and are degenerate 
in the primary CMB measurement. In Fig. |7| we forecast the accuracy for ACT to measure the kinetic SZ C; assuming the 
reionization scenario B. Al — 2/10 is adopted. At I < 2500, primary CMB dominates and prohibits the extraction of the 
kinetic SZ signal. At I > 10 4 , the instrumental noise dominates. At I ~ 6000, ACT will reach its highest sensitivity with a 
3% accuracy in the kinetic SZ Ci measurement. Since the kinetic SZ power spectrum is almost flat at I > 3000, one can bin 
the data points (~ 10 points, fig. |7J at 3000 < I < 8000 and do not lose any crucial information. Binning in this way, the 
ACT-measured Ci accuracy will reach better than 1%. 

The accuracy of inferred reionization history depends on both the precision of KSZ power spectrum measurement and theo- 
retical predictions. Given perfect theoretical predictions, for the single reionization scenario, the uncertainty in the reionization 
redshift Az ro i on /.z rc i on ~ AC;/Ci/0.34 ~ 3%. The two CMB-degenerate reionization scenario B & C can be distinguished with 
> 10a confidence 2 . But uncertainties in the KSZ signal extraction and theoretical prediction may degrade the power of ACT 
to measure the kinetic SZ effect. 

One is the lensing effect. Since lensing changes CMB photon positions but not CMB temperature, it results in a smoothing 
of the CMB power spectrum. Because the primary CMB power spectrum drops quickly after I ~ 2000, lensing amplifies it by 
large factors at small angular scales (Fig. . In our estimation, we have assumed that this lensed CMB power spectrum can 
be predicted and subtracted. Since the lensing effect at these scales involves non-linear dark matter clustering, the precise 
prediction of its power spectrum could be difficult. If we take this possibility, the residual CMB (lensed subtracted unlensed) 
would prohibit the measurement of KSZ at I < 4000 and increase the error of KSZ Ci measurement at I £ [4000, 6000] to 
be larger than 10%. But its effect at I > 6000 is small. Binning the ~ 4 data points at 7000 < I < 9000 would still give a 

2 We assume different bins are uncorrelated. This assumption is exact for a whole sky survey. For a survey with limited sky coverage, 
this assumption only roughly holds in our case. For two successive data points of Fig. [7|at 3000 < I < 8000, the corresponding spacial 
separation at z ~ 10 is ~ 1 Mpc/h, which is about the same as the density correlation length at z = 10. So, one can safely neglect cross 
correlations in different bins. 
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better than 2% accuracy. Furthermore, the future lensing survey such as the lCFHT legacy survey! <l2003l) will measure the 
dark matter distribution to 1% accuracy and allow a precise prediction of the lensed CMB power spectrum. So, the lensing 
effect of the primary CMB could be well handled. Lensing also changes the KSZ power spectrum. But since the KSZ power 
spectrum is almost flat, this lensing effect is negligible. 

In our discussion, we have assumed t he thermal SZ effect to be negligible at ~ 225 Ghz. This is true in the non relativistic 
limit. The relativistic correction (see e.g. iDoleov et all I^OOlh - ) p roduces a temperatur e decrement at this band and its effect 
on ACT measurement of cluster SZ effect has been discussed llAghanim et alJ 12002'). This relativistic correction is of the 
order k B T e /(m e c 2 ) ~ T e /500keV. Though it is non-trivial for the measurement of high temperature cluster SZ effect, it is 
negligible for a blank sky SZ power spectru m measurement. The mean elec tron temperature at present due to gravitational 
heating is expected to be less than 0.5 keV JPerJll99St Izhang fc PenlbOQll) . COBE/FIRAS puts an upper limit 1.5 x 10~ 5 
jFixsen et al. Hl99rJ) on the mean 'y' parameter and basically excludes dramatic feedback with T > 1 keV. One then expects 
the mean electron temperature at present to be < 1 keV. Sin ce the blank sky the rmal SZ power spectrum at I > 4000 is mainly 
contributed by these low temperature IGM gas at z > 0. 5 fehang fc Penll200lh . we expect that the relativistic correction is 
much less than 1% and is thus negligible. 

We h ave not considered the con tamination of unresolved IR and radio sources, whose effects are estimated to be non- 
negligible ijwhite fc MaiumdaJl200^ . At small angular scales, these sources are present as Poisson noise in the KSZ mea- 
surement. From the large I power spectrum behavior, the power spectrum of the sum of instrumental noise and this Poisson 
noise can in principle be measured and subtracted. Multi-frequency observations also help to clean them out. 

We estimate the dependence of the KSZ power spectrum on cosmological parameters. One can easily work out that 
A|(fc,z) « p 2 H 2 {z)al~ 6 (4 in the linear regime and 6 in the highly nonlinear regime where stable clustering holds). Here, 
13 = adD/da/D ~ Q™(z). In our interested 3000 < I < 10000, Ci contributions mainly come from z > 1, where j3 ~ 1, 
H(z) oc SIq and the comoving distance x is roughly oc fl^ 1 with < 7 < 1/2. Here, Qo = Q. m {z — 0). Applying the above 
relation to Eq. [3] we find that C ; oc {Q, b h) 2 Q} /2+l(nB ~ 1) <ri~ 6 oc (ft b h) 2 fl~ 1/2 at 6 . As a reminder, riB is the power index of 
A 2 3 (k, z) at k = l/x and < nj < 1 for I > 10 3 and z > 1. So there is a degeneracy between jz rc i on , Q b h and <rs and the 
recovered z rc ion oc (Q b h) 6 tjg 15 . Thus, current uncertainties in cosmological parameters degrade the accuracy of the extracted 
Zrcion significantly. But this degenerancy is not that bad as appears since it is the matter density power spectrum PjvL(fc) at 
Mpc/h scale, instead of erg, that directly determines the kinetic SZ effect at I ^several thousand s through Ci oc Pfj L . Pnl {k*) 
at Mpc/h scale will be measured by CFHT legacy survey with 1% accuracy (see Pen et alJ l|2003h for the discussion of relevant 
length scales in the recovered 3D power spectrum). Future experiments such as Planck, S NAP and CFHT would signific antly 
reduce the above uncertainties in cosmological parameters and density power spectrum jTegmark fc Zaldarriagall2002l) and 
allow a precise recovery of the reionization history from the kinetic SZ effect. 

At the presence of various noise sources in KSZ measurement discussed above, ACT would still be able to measure 
Ci with ~ 2% accuracy. Given a sufficient understanding of cosmic baryonic matter density and the matter density power 
spectrum, which will be measured precisely by future experiments, the reionization redshift could be constrained with several 
percent accuracy, some multiple reionization scenarios degenerated in the primary CMB temperature and T-E polarization 
measurement could be distinguished by ~ 10<r confidence. 

By far we only discussed the uniformly reionized universe. Patchy reionization introduces extra power to the KSZ effect 
and may degrade the power of KSZ to constrain the reionization history. But future observations may be able to separate 
it by its characteristic Ci behavior. Furthermore, by cross correlating the KSZ effect with the cosmic 21 cm background, 
one can further recover the redshift information of the KSZ effect at high redshifts and extract the contribution from th e 
patchy reionization. S i nce in the literature , the significance of patchy reionization is under debate fe.g. lGnedin fcJaffel J200lt) : 
IValageas et alJ l)200lh : lsantos et alJ ^2003^ 1. we postpone its further discussion. 



8 SUMMARY 

The dependence of the kinetic SZ effect on the reionization history allows a detailed probe of the reionization process and 
the feasibility to break the degeneracy posed in the primary CMB measurement. This relies on the precision of the kinetic SZ 
effect measurement, which will be realized by future CMB experiments such as ACT and precision modeling of the kinetic SZ 
effect. We have built an analytical model for the kinetic SZ effect and tested against our high resolution hydro simulations 
through the 3D power spectrum A%(k, z) of the momentum curl part, which completely determines the kinetic SZ effect 
angular power spectrum. The grid version of our model, which captures various limitations of simulations such as the finite 
box size effect and thus reproduces simulation results at all available simulations scales and redshifts with better than several 
percent accuracy. The continuous version, which represents an ideal simulation with infinite box size and resolution and is 
free of simulation limitations, should describe the real A 2 B to a high accuracy. Comparing between simulations and analytical 
predictions, we found that 

• The simulated A 2 B loses power at large scales due to the cutoff of the density fluctuation beyond the simulation box. 
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This effect causes the simulated A| to be steeper than the real one at large k. The smaller the simulation box, the severe the 
problem is. 

• The simulated A% loses power at small scales due to the finite resolution. 

• Correspondingly, the simulated kinetic SZ power spectrum loses power at both large and small scales. In order for 
simulation to directly produce reliable Ci at fOOO < I < 1 0000 without reference to our models, a hydro simulation with box 
size > 400 Mpc/h and resolution > 4096 3 is required. 

• The effect of the non-linearity on A 2 B has the same tendency as its effect on the density power spectrum. Its effect can 
be modeled by substituting both linear density power spectra in the Vishniac integral (Eq. |?J & |SJ with the corresponding 
non- linear ones (Eq.|Hl&|7J. For WMAP cosmology, the non-linearity increases the density power spectrum and A B . 

• Our analytical model in the grid version captures all main features of simulated A% and reproduces the simulation results 
with high accuracy. Its continuous version is an accurate description of the real A B and therefore the KSZ effect over a wider 
parameter range in power spectrum slope than the past models proposed in the literature. 

• For the WMAP cosmology and the single reionization scenario, the kinetic SZ power spectrum can be approximated with 
better than several percent accuracy as l 2 Ci/(2n) ~ 0.91 x 10" 12 [(1 + 2 rcion )/10] - 34 (;/50 00) - 23 " ' 015(Zrt!i °"" 9) for 6 < z Ieion < 
20 and 3000 < I < 9000. This strong reionization redshift dependence will allow a precision measurement of the reionization 
history. For the WMAP favored z Te i on ~ 16.5, the kinetic SZ power spectrum has a broad peak extending from I ~ 3000-10 5 
with an amplitude ~ 2.7/ik. Among cosmological parameters, ag and Qb are the two that C; is most sensitive to. C; scales as 

(n b h) 2 at 6 - 

• ACT will measure the kinetic SZ effect to 1% accuracy. Given precise measurements of cosmic baryonic density Qbh 
and matter density power spectrum, for the single reionization scenario, the reionization redshift z Ic i on can be constrained 
with several percent accuracy. The kinetic SZ effect can further distinguish some more complicated reionization scenarios and 
breaks the degeneracy met in the primary CMB measurement with ~ 10cr confidence. 
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